function h = plot_u(grid,base,U,varargin)
% h = plot_u(grid,base,U)
% h = plot_u(grid,base,U,k)
% h = plot_u(grid,base,U,k,h)
%
% Plot a scalar U (represented as a discontinuous field, i.e. as a two
% dimensional array). On each element, the solution is resampled on a
% local mesh, which is then plotted with my_pdesurf. If k is absent,
% the degree of the basis is used.

 % 0) preliminary checks
 if(nargin>=4)
   k = varargin{1};
 else
   k = base.k;
 end
 if(nargin>=5)
   h = varargin{2};
   figure(h);
 else
   h = figure;
 end
 
 % 1) resample the solution
 [Xi,Ti,Ui] = el_resample(grid,base,U,k);

 % 2) open the figure and plot the mesh
 my_pdemesh(h,grid,'k');
 hold on

 % 3) plot the field U
 switch grid.d
  case 1
   for ie = 1:grid.ne
     my_pdesurf(h,Xi(:,:,ie),Ti,Ui(:,ie));
   end
  case 2
   np = size(Xi,2);
   nt = size(Ti,2);
   XXi = zeros(2,np*grid.ne);
   TTi = zeros(4,nt*grid.ne);
   UUi = zeros(np*grid.ne,1);
   for ie = 1:grid.ne
     XXi( : ,(ie-1)*np+1:ie*np) = Xi(:,:,ie);
     TTi(1:3,(ie-1)*nt+1:ie*nt) = Ti + (ie-1)*np;
     UUi(    (ie-1)*np+1:ie*np) = Ui(:,ie);
   end
   TTi(4,:) = 1;
   my_pdesurf(h,XXi,TTi,UUi);
  otherwise
   error('This functions supports only 1D and 2D.')
 endswitch
 axis equal

 hold off
 
return

